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Abstract 

We study kinetic model of reversible enzyme reaction and compare 
two techniques for analytic approximate solutions of the model. Analytic 
approximate solutions of non-linear reaction equations for reversible en- 
zyme reaction are calculated by Homotopy perturbation method (HPM) 
• and Simple iteration method (SIM). The results of the approximations 

. y are similar. The Matlab program are included in appendixes. 

Keywords Enzyme Kinetics, Homotopy perturbation method, iteration method, 
Q" 1 Michaelis-Menten Kinetics, Quasi-steady state approximation. 

1 Introduction 

The variety of chemical reactions in a living organism is carried out by enzymes. 
f-L It seems to be the rate of chemical reactions (both forward and backward) is 

accelerated by enzymes. They are very vital because many chemical reactions 
qq are occurred without the activity of enzymes. They are bound by an enzyme's 

active site and they become a product after a series of stages. These stages are 
known as the enzymatic mechanism. There are two types of mechanism, single 
substrate and multiple substrate mechanisms [T7J [TTJ [251 HI]- An important 
branch of enzymology is enzyme kinetics and these kinetics are used to study 
of the rate of chemical reactions. Some differential equations are involved to 
characterize the enzyme kinetics with some principle of chemical kinetics jlOj . 
To comprehend the role of enzyme kinetics, the researcher has to study the rates 
of reactions, the temporal behaviors of the various reactants and the conditions 
which influence the enzyme kinetics. Introduction with a mathematical bent is 
given in the books by Rubinow [IB], Murray [13] and Segel [T§] . 

The single enzyme reaction is one of the most powerful area that many 
attempts have been shown to find what is the behavior of this reaction based on 
their parameters, in the following paragraphs we try to show some of previous 
works regarding this reaction. Simply this enzyme reaction is defined as follow, 

E + S^ES^P + E (1) 
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Where [E] is Enzyme concentration, the substrate concentration is defined by 
[S], [ES] is Enzyme-substrate complex and [P] is product. Also, ki,k 2 and k3 
are represented the reaction rate constants. By using the idea of mass action 
simply we can show the reaction equation ([T]) in terms of system of non-linear 
ordinary differential equations [21] . 

There is a variety of possible simplifications for system to describe analytic 
approximate solution of it. One of the most common attempt to simplify this 
system is qassi-steady state approximation (QSSA). The quasi-steady state as- 
sumptions are occurred as a fundamental assumption for enzyme kinetics and 
the history of this subject began 80 years ago. It has a grate role to analysis 
the equations of enzyme kinetics [TU]. And also, it has many applications on 
the field of experiment results. There are two laws of quasi-steady sates. First 
law is that the enzyme-substrate complex concentration just occurs in terms of 
constant rate. The second one is that the concentration of enzyme-substrate 
complex is changed in the small rate. 

Another simplification is returned to Michaelis and Menten in 1913 who pointed 
out the enzyme reaction ( 1 1 should be k 2 >• k 3 , therefore = jj^ . It means, 

there is an equilibrium between [E] , [S] and [ES] to produce [P] and [E] . 
Briggs and Haldane proposed that in 1925 Michaelis and Menten assumption 
is not always applied and they said " it should be replaced by the assumption 
that [ES] is present not necessarily at equilibrium but in a steady state under 
condition [So] ^ [Eq] " . Also, according to their assumption the rate of enzyme- 
substrate complex [ES] and its rate of consumption during the reaction. It 
means, the concentrations [ES] is occurred as a steady state. This is known 
as the steady state assumption (SSA) or sometime is called quasi-steady state 
approximation (QSSA), or seduo-steady sate approximation [7]. 
The first description of QSS was defined by Briggs and Haldane in 1925 pp. 
They describe the simplest enzyme reaction in equation (JlJ and they pointed 
out the total concentration of enzyme [E], where [E] to t = [E] + [ES] is a tiny 
value in comparison with the concentration of substrate [S]. Also, they have 
shown the term of d ^f^ is negligible in compared to and As a result, 
they found the Michaelis-Menten equation which is a differential equation used 
to describe the rate of enzymatic reactions. The classical Michaelis-Menten 
equation is defied as follows, /ci = (k 2 + ks)[ES], or 

^ = ^™ = S?1 

where Um = fc2 ^ fca is the Michaelis-Menten constant (for more details see [5] ). 
There is many incorrect comment in some papers which were published about 
QSS [5]. Two of them are: 

1. Fast reactions include the concentration of enzymes. As a result, the 
concentration of substrate or stable components are slower than relax. 
This result basically is not true for some QSS systems. For instance, 
the majority of reactions in the Michaelis-Menten consist of enzyme with 
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product or substrate. It is difficult to find any fast reaction about enzyme 
without product or substrate. 



2. Sometimes in QSS systems the time derivatives equate to zero, hence the 
concentration of intermediates occur as a constant. Generally speaking, 
this is obviously not correct. Because of this, we can equate the kinetic 
reactions for some time derivatives to zero and these derivatives usually 
are not equal to zero when we apply the QSS for the system. To show this 
fact we can easily differentiate the equation ^ , then we can find [ES] for 
QSS. It is a constant when [E] 3> Km and this condition is different from 
the Briggs-Haldane condition [E] + [ES] < [S]. 

Shnell and Mendoze (1997) based on stander quasi-steady state (sQSSA) have 
shown a solution for total time evolution for the substrate [S] for equation 0, 
that is 

[S'](t) = W([S ]exp(-kt+[S ])), (3) 

where [S'j = jp-, k = v m axl Km and W is satisfied the transcendental equation 
W(x)exp(W(xf) = x. 

After that Laider (1955) pointed out in [T7] that the actual validity of (QSSA) 
and he derived that the substrate [S] greatly exceed of the enzyme [E] . It means, 

^1 « l (4) 
[So] 1 ' 

Recently, Segel (1988) and segel and Slemrod (1989) in [17] emphasized that 
the more general validity condition for (QSSA) is defined as, 

[E ° ] « 1 (5) 



(K M + [S ]) 



As a result, these conditions have a great impact to find many attempts to ap- 
proach the best approximate solution to equation and to explain different 
cases for this simple chemical reactions. Another attempt, Schenell and Maini 
(2000) proposed that the reverse quasi-steady state approximation (rQSSA) as- 
sumption, they assumed that there is a stationary value in terms of [ES] for 
equation ([T]) at high enzyme concentration |17j . For this reason, this assump- 
tion has become a basic in enzyme kinetics and it has applied to find analytic 
approximate solution and explain the parameters of many chemical reactions. It 
seems that one of the powerful application of (QSSA) in the field of biochemical 
system is used to reduce the parameters and variables in a complex system in 
to simple one. This kind of simplification is useful in biological systems such as 
genetic regulation process and metabolic process. Another important applica- 
tion for (QSSA) in terms of data information from experimental is realize that 
the (QSSA) and Micheales-Menten equation have an important role to describe 
of enzyme kinetics whenever the concentration of substrate largely exceeds that 
of concentration of enzyme. 
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More recently, Tzafriri showed that the total quasi-steady state approximation 
(tQSSA) for equation ([!]) has uniformly validated by self-consistency for initial 
transient [20] . 

After that, Morten has used the (tQSSA) for some complex reactions such as 
commutative system and the Goldbeter-Koshland switch reactions [2]. He has 
shown that the (tQSSA) provides the best approximate solution to the systems 
for a wide range of parameters and he emphasized that the (tQSSA) is much 
more better than (sQSSA) to estimate parameters and to find real value of 
them. The techniques are used by many researchers to reduce the number of 
parameters and variables in the field of biochemical reactions, but it could be 
said that sometimes the behaviors of a new system are quite different compared 
to original one, especially to estimate stability of chemical reactions, this is one 
of the negative point of the reducing the parameters and variables [3] . 
Non linear models with multiple timescales are described by Quasi-steady state 
(QSS) and Quasi-equilibrium (QE) approximation. They are normally distinct 
physically and dynamically [IB] . Gorban and Shahzad in [F] have described 
complex chemical mechanisms under two assumptions which are the Quasi- 
steady state hypothesis (QSS) and the quasi-equilibrium hypothesis (QE). They 
have proved the generalization of mass action laws in terms of basic relations 
among kinetic factors. Briggs and Haldane proposed the (QSS) for equation ([!]) 
and they showed that the total concentration of enzyme [E] tot is much lower 
than the total concentration of substrate [S]tot, where [E] to t — [E] + [ES] and 
[S]tot — [S] + [ES], therefor the low concentration, fast species of equation ([I]) 
is the concentration of [ES]. In addition, the concentration [S] has a grate 
role to drive [ES]. As a result, the reaction ([I]) can be simplified for a unique 

irreversible reaction S R ■ s _i_^ tot ' ) p > means that = — ^1 = k 3 [ES]Qss- 
The QSS for complex concentration [ES] comes from the equation, 

kx[S]([E] tot - [ES] tot ) = (k 2 + k 3 )[ES] QSS 

So we can rearrange the above equation in terms of [ES]qss as follows, 

k 2 [E] t ot[S] 



R{[S],[E] tot ) 



k m + [S] 



The QE approximation has been applied for the system equation ([I]) which the 
first reaction is fast, reversible reaction. The simplification of equation ([I]) in 
this case based on a pooling of species. The fast reversible reaction lead to 
reserve two pools [E] to t and [S] to t- Furthermore, the two reactions of the sys- 
tem have a grate role to consume [-E]t t- Whereas, the second reaction in the 
mechanism consume the pool [S 1 ]^. It is called the slow variable of the system. 

The single step approximation is defined by Stot^ ^ ''"''-P, or equivalently 
d\P\ _ _d[^a 1 _ k 3 [ES]QE- The positive solution of the quadratic equation, 
k 1 ([S] tot -[ES]QE)([E] tot -[ES] QE ) = k 2 [ES] QE is the value of QE of the com- 
plex [ES]. For this, w e can find that, R([S] to t[E]tot) = 2k 3 [E] tot [S] t ot([E] tot + 
[S]tot + If)" 1 , (1 + V 1 ~ 4[E]tot{S]tot/([E]tot + [S] tot + feAi) 2 )- 1 , when the 
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total concentration of enzyme is much small than the total concentration of 
substance. It means that [E] tot <§C [5]t t. So, we obtain the basic equation of 
Michaelis and Menten, R{[S] toU \E] tot ) a fc 3 { f^} tot ■ 

The purpose of this work is to derive asymptotic approximate expressions for 
the substrate, product, enzyme and enzyme-substrate concentrations for equa- 
tion ^ by using Homotopy perturbation method (HPM) and Simple Iteration 
method (SIM) and to point out the similarities and differences between of them 
for all values of dimensionless reaction diffusion parameters e, A, a and k. An- 
other aim for this project is to find out the appropriate iteration in (SIM) which 
is fitted to our solution. 



2 Mathematical Formulation 

The equation of of Michaelis and Meten has been applied. It had been 
succeed by Kuhn in 1924 [12] to several cases of enzyme kinetics. The model 
of biochemical reaction was developed by Briggs and Haldane in 1925 [3T] . The 
model of an enzyme action considers a reaction that includes a substrate [S] 
binds an enzyme [E] reversibly to a substrate- enzyme [ES]. The substrate- 
enzyme leads reversibly to product [P] and enzyme [E\. This mechanism is 
often written as 

E + S^ES^P + E (6) 

This mechanism shows the binding of substrate [S] and release of product [P] 
where the free enzyme is [E] and the enzyme-substrate complex is [ES], also 
ki,k2, &3 and k 4 denote the rate of reactions. It is clear that from equation (|6| 
substrate binding and product are reversible. The concentration of the reactants 
in the equation ^ is denoted by lower case letters 

e=[E],s = [S],c=[ES],p=[P] (7) 

The time of evolution of equations ^ and ^ is found by the law of mass action 
to obtain the set of system of following non-linear of reaction equations. 

ds . . 

— = -fcxes + k 2 c (8) 

de 

— = — kies + (fc 2 + k 3 )c — k^pe (9) 



dt 
dc 

n 



k\es — (}t2 + k^)c + k^pe (10) 
-5 = fc 3 c- k 4 pe (11) 



when the initial conditions at t = are given by 

e(0) = eo, s(0) = s , c(0) - 0,p(0) = (12) 
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Adding equations Q and (10), and using initial conditions (12), we obtain 

e + c = e (13) 

Also, adding equations Q, (10) and (11), and using initial conditions (12), we 
get 

s + c + p = s (14) 

By using equation (13) and equation (14), the system of ordinary differential 
equations ([8|-(11) reduce to only two variables s and c as follows, 



— = -feie s + (fcis + k 2 )c 
at 



dc 
dt 



— kie s + (k 1 s + k 2 + k 3 )c + fc 4 (e — c)(s — s — c) 



(15) 
(16) 

UL 

ith initial conditions s(0) = Sq, c(0) = 0. By introducing the following param- 

u\ „r+\ ^(+\ (n\ 

,u(t) = 

(17) 



with 
ctcrs 



r = 



•so 



eo 



A = 



k 2 , k 2 + k 3 e fc 4 

— ; ,£= — , a = — , m = A + ae + a 

so 



,fc = 



We use the 



fciS ' kiS 

dimensionless technique to reduce the number of parameters for the 
nna+irvriQ /TiKh QnH lupii an d the initial conditions (12) and it can be 



system of equations (15) and (16) 0,1. w . .. 
represented in dimensionless form as follows: 

du 



— = — eu + e(u + k — A)v 

(XT 



— = it — hi + k)v + all — v)(l — u — ev) 
dr 

dw o 

— — = cm — + ev +1) 

AT 

i(0) = l,v(0) = 0,w(0) = 



(18) 

(19) 

(20) 
(21) 

In this paper we estimate the analytic approximate solution for system of non- 
linear ODE equations ((TsJ-Q, by using the methods of (HPM) and (SIM). 



3 Analytical Approximate Solution by Homo- 
topy perturbation method 

Basic idea of Homotopy-Perturbation method (HPM) is defined in this section 
and then it applies to find the approximate solution of the problem equations 



(18)- (21). It is considered from the following function 

A{x) - f(r) =0,re(! 



(22) 
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with the boundary conditions 



B(x 



dx, 
dn' 



0, r G r 



(23) 



where A, B, f[r) and are general differential operator, boundary operator, a 
known analytic function and the boundary of the domain SI, respectively [8]. 
The function A consists of linear part L and non-linear part N. So the equation 
( 22 ) can be written as, 

L{x) + N(x) - f(r) = (24) 
The Homotopy function is defined by z(r,q) : SI x [0, 1] — > R, which satisfies 

H(z,q) = (1 - q)[L(z) - L(x )] + q(A(z) - f(r)) = 0,p € [0, 1], r € SI, (25) 

or, 

H(z, q) = L(z) - L(x ) + qL(x ) + q[N(x) - f(r)} (26) 

where q € [0, 1] is an embedding parameter. Whereas Xq is an initial approxima- 
tion of equation (22), which is satisfies equation (23). Basically, from equation 
( 25 ) and equation ( 26 ) we can obtain 

H(z,0)=L{z)-L(x o ) = 0, (27) 

H(z,l)=A(z)-f(r) = (28) 

Changing z(r, q) from xq to x(r) depends on the values of q from zero to unity. It 
is called deformation in field of Topology, whereas L(z) — L(xo) and A(z) — f(r) 
are called Homotopy. We use q as a small parameter at the first time and we 



defied the solution of equation (251 and equation \2Q\ as a power series in q: 

(29) 



z = z + qzi + q z 2 + 



Let q = 1 to get the approximate solution of equation ( 22 ) 

zo 



lim z 

1 



+ Zl + Z 2 + 



(30) 



Thus, HPM includes the combination of perturbation method and the Homotopy 
method. The equations ( 18 )-([20| can be solved analytically in a simple and 
closed form by using HPM in (Ref appendix A). So, the approximate solutions 
of the system of non-linear differential equations ( 18 1 and ( 19 ) become, 

ab ae\ _^ abc—aae + aac 



u{t) =2e- £T + 



as 



as 

1 

C — £ C 

cbe — cae 



te- £T 

(-e-c)t i 



c 2 (e — c) 
aae - 

1 + 



abc 

c(e — c) 2 



aca 



act 
ce 
aa 
ce 



c(e — c) 2 
b 



-2er 



£ — C 

b 



-as 



cbe 



c 2 (e 



(31) 
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v(t) 



be — as + ca _ c 
c(e - c) 
ba 



+ 



c — e c(c — e) 
aeb — cb 2 — cba 



-2et 



ce{e — c) 
b 2 

(e-c)(c-2e) 



(-e-c 

' v £ — C 

cfe 2 — bae + c6a 

ce(e — c) 



(c- e)(c- 2e) 
6 



Q 
C 



(32) 



c(e — c) 



The analytic expressions of the substrate u(t) and enzyme substrate v(t) con- 
centrations can be represented in equations (pll) and (32 1. The dimensionless 



concentration of enzyme E can be obtained from the equations (13) and (17 1 
as follows, 

E(t) = ^l = l- v(t) (33) 
eo 

The dimensionless concentration of the product w is obtained either by equation 



(20) as follows, 



w (t) = J (a(u(t) - u(t)v(t) - ev 2 (t) - 1) + mv(t))dt 



(34) 



or, we can use equation (17) and equation (14)to find the concentration of the 
product w as follow 



w(t) = 



1 — u(t) — £v{t) 



(35) 



The simple analytic approximate solution form of the concentrations of enzyme 
E(t) and product w{t) for all values of parameters e, A, a and k are represented 
in equations (33l-([35|. We can easily realize the behavior of the solution based 
on this method (HPM) from the Figures [l6|20j 



4 Simple Iteration Method 

In this section we use a simple technique that is new to find the analytic ap- 



proximate solution for the system of equations ( 18 ) and ( 19 ) and we introduce 



this method by rewrite the equations ( 18 ) and ( 19 1 as follows 



du 

7k 



-eu + e(k — X)v + euv 



dv 
7h 



= (1 — a)u — (k + a + ae)v — (1 — a)uv + aev 2 + a 



(36) 



(37) 
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Let a = e(k — A), b = (1 — a) and c = k + a + ae then the equations (36) and 
(37 (can be written as: 



"n+l 



A 



Un+l 
V n +1 



G{u ni v n ) 



for n = 0,1,2,... where, G{u n v n ) 



eU n V r , 



aev„ 



(38) 



is a non-linear 



part of system ( 38|) , and 

-e a 



A = 



b -t 



is a matrix for the linear part of equation ( 38 ) . To eval- 



uate an approximate solution of equation ( 38 ) with initial conditions equation 



(21), we introduce the following steps to approach the approximate solution: 



Step 1. For n — 0, Uq(t) = 1, Vq(t) = and if possible suppose that a — > 



just in this step. It means that we assume non-linear part of equation (38) 



approach to zero, by putting in equation ( 38 ) , we obtain the following system 



.4 



Ml 



(39) 



We can solve the system of ordinary differential equations 



So, the solution of equation (39) with initial conditions (21) 



(39} analytically [2]. 
is 



ui(r) 
vi(t) 



d 2 e PlT 
d ie piT 



d 3 e P2T 
die P2T 



(40) 



where pi and p 2 are eigenvalues of the matrix A, and d\ 



(P2+e) 



«1e 



-Pi) 



and d?, = 



(Pi+e) 

(P1-P2) ' 



_ (Pl+£)(P2+g) 
0(P2-Pl)_ 



, d 2 



We substitute u\ and v\ in equations (33) and 



(35), then we obtain Ei and Wi, respectively. The behavior of concentrations 
of equation (40) are described in Figures l][5 (see Appendix C). 



Step 2. For n = 1 and substituting equation (40) in equation (38), we obtain 
the following system of non-linear ODE, 



= A 



"2 



G(ui,ui) 



(41) 



v 2 J \ v 2 

We can easily use methods for solving the system of non-linear differential equa- 



conditions equation (21) is 



tions to solve the system ( 41 ) [2] . The solution of above system with initial 



u 2 (t) = ac 3 e PlT + ac 4 e p * T + d 22 e 2piT + d 23 e^+ p ^ T + d 24 e 2p2T + d 25 (42) 

v 2 (t) - c 3 (pi + e)e p ^ + Ci (p 2 + e)e p * T + d 26 e 2piT + d 27 e^ +p ^ T + d 2S e 2p2T + d 29 

(43) 

Where d 22 ,...,d 2 g and 03,04 are constants. We substitute u 2 and v 2 in equa- 



tions ( 33 ) and (35 1, then we obtain E 2 and w 2 , respectively. The behavior of 
concentrations in this step are described in Figures 6||To| (See Appendix D) . 
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Step 3. For n = 2 and substituting equations (42) and (43) in equation (38) 
we get the following system of non-linear ODE, 



= A 



«3 



G(u 2 ,v 2 ) 



(44) 



We can also use methods for solving the system of non-linear differential equa- 



tions to solve the system (44) [2J. The solution of above system with initial 
conditions (21 ) is 



u 3 (t) =ac 5 e PlT + ac 6 e P2T + d 9Q e 2piT + d 91 e ipi+P2)r + d 92 e 3piT 

+ d 93 e {2pi+P2)T + d 9i e {pi+2p2)T + d 95 re PlT + d 9e e PlT 

+ d 97 e 2p2T + d 98 e 3p2T + d 99 e P2T + d W0 re P2T + d wl e 4p ' T 

+ d 102 e (3pi+P2)T + d 103 e (2pi+2p2)T + d 104 e (pi+3p2)T + d 105 e ip2T + d 106 

(45) 



v 3 (t) =c 5 h ie pir + c 6 h 2 e P2r + d 107 e 2piT + d 



(pi+P2)t 



d\\ 2 Te 



Pit 



lose 

d W9 e 3piT + d lw e (2pi+P2)T + d lll e {pi+2p2)T 
dn 3 e PlT + d lu e 2p2T + d 115 e 3p2T + d 11& e P2T 
d 117 re P2T + d 118 e 4pir + dn 9 e (3pi+P2)T + d 120 e (2pi+2p2) 



(46) 



+ d 121 e^ +3p2 ^ T + d 122 e 4p2T + d 123 
Where d 99 , ...,di 23 and c§,c§ are constants. We substitute u 3 and ^3 in equa- 



tions ( 33 ) and ( 35 1 , then we obtain E 3 and w 3 , respectively. The behavior of 
concentrations of this iteration are described in Figures TTp5 (see Appendix E). 
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5 Asymptotic Analysis 



Asymptotic analysis was suggested by Kruskal (1963) for the field of chemical 
reaction [6]. He defined asymptotology as " the art of describing the behavior 
of a specified solution (or family of solutions) of a system in a limiting case" . 
The initial concentration of enzyme E and concentration of substrate S led to 
a great effect on the concentration of free enzyme E and substrate S. So, it can 
be identified the following three different conditions based on initial ratio y|4 

El- 

a) If the initial concentration of enzyme [Eq] is much more great than the 
initial concentration of substrate [So]. It means, j^j ^> 1. Also, according 
to S. Schenell and P. K Maini in [TS] emphasis that the initial concentration of 
enzyme is greatly exceeded of the concentration of substrate, that is [E ] 3> [So]. 
So, from equation (17) we can point out e > 1. In this case, the part of enzyme 
concentration which is bound to the concentration of substrate is small. It 
means, there is a free rate of enzyme, this rate based on the availability of 
the substrate and it is increased whenever the concentrations of substrate is 
increased or by adding the extra of substrate to the chemical reaction. 

b) If the initial concentration of substrate [So] is much more grater than the 
initial concentration of enzyme [E ]. It means, j^j <C 1. So, from equation 
(17) we obtain e < 1. In this case, there is a small part of of substrate which 
is bound to the enzyme and a part of it is free. While, enzyme molecules in 
this case usually bound to substrate molecules. It means, there is a small rate 
of enzyme which is free. The availability of enzyme in this case depend on this 
rate and increase when the rate of enzyme is increased or by adding some extra 
enzyme to the chemical reaction. 

c) If the initial concentration of e nzyme [Eq] and substrate are equal. It means, 
|^| = 1, so from equation (17) we get e = 1. In this case, there is no any 
free molecules of enzyme and substrate. It means, all substrate molecules are 
occupied by the molecules of enzyme and all enzyme molecules are also limited 
by the molecules of substrate. 

Furthermore, in terms of the rate constant of reactions k^ and k\ from equation 
(17), we can define the following conditions: 

d. If foj ^> ki then a > 1 . 

e. If k\ <C k\ then a < 1 . 

f. If &4 ~ fci then a ~ 1 . 

In addition, according to definition of A and k from equation (17), we obtain 
just one condition that is A < k because k^ is always positive value. Generally 
speaking, we can describe the behavior of analytic approximate solutions of u, v, 
E and w in terms of different cases of the value of parameters. As result, we can 
easily combine the condition(a)-condition(f ) to find all possibilities regarding the 
parameters e and a. This combination consists of five different cases. These 
cases not depend on the parameters A and k because it is clear that A < k from 
equation (17). Thus, we can define the following five basic cases in this paper: 
Case 1. The value of e ~ 1 and a ~ 1. 
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Case 2. The value of e > 1 and a > 1. 

Case 3. The value of e < 1 and a < 1. 

Case 4. The value of e > 1 and a < 1. 

Case 5. The value of e < 1 and a > 1. 
With A < fc for each case. As a result, we can apply the various values of 
parameters to both methods (SIM) and (HPM). It seems that these cases are 
very helpful to make distinguish between the methods which are used for our 
problem and it could be the best way to choose the appropriate iteration of 
(SIM), in comparison with (HPM). Another positive point about the asymp- 
totic is that they have an important role to combine information between math- 
ematical model and experimental data. Therefore, the cases tell us the simple 
approaches which are used in this work are applicable and easy to describe the 
chemical reaction behavior of equation Q in terms of our five cases. As a re- 
sult, we apply the above cases to show the result of our solution to point out 
the similarities and difference between the (HPM) and (SIM). 



6 Results and discussion 

Figures [T]|20] in this section show the analytic approximate solution of substrate 
u, Enzyme E, enzyme-substrate complex v and product w. Our results based 
on two methods which are mentioned in previous sections. Simple iteration 
method (SIM) is used to describe the approximate solution in Figures [TflT5l and 



Homotopy Perturbation method (HPM) is also used to describe the approxi- 
mation solution in Figures [l6p0| for the same problem. We can normally point 
out that the figures change in terms of values of dimensionless parameters e, a, 
A and k. Each figure in this work corresponds to one case from case 1-case 5 
in previous section. It is more easy to estimate the similarities and differences 
between the two methods which have been used in our work. Furthermore, it is 
clear that changing the value of parameters leads to obtain different value of the 
substances. The figures in this section are divided in to four groups of figures, 
the first three groups are related to three iterations of (SIM) and the last group 
of figures refer to the (HPM). 
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Figures [T](5] in these profiles of the normalized concentrations of the substrate 
u, enzyme-substrate complex v, enzyme E and product w correspond the case 
1-case 5, respectively. The equations of step 1 are applied to plot the Figures 
(see Appendix C). 




01 234 5 6789 10 
Dimensionless Time{t) 



Figure 1: e = 1, a = 1, A = 0.4 and k = 1.3 
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Figure 2: e = 1.6, a = 1.3, A = 0.9 and k = 1.7 
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Figure 3: e = 0.8, a = 0.2, A = 0.6 and k = 1.1 




01 23456789 10 



Dimensionless Time(t) 

Figure 4: e = 1.3, a = 0.3, A = 0.9 and k = 1.2 
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Figure 5: e = 0.6, a = 1.2, A = 1.2 and fc = 1.7 

Figures [6] jl0| in these profiles of the normalized concentrations of the sub- 
strate u, enzyme-substrate complex v, enzyme E and product w correspond 
the case 1-case 5, respectively. The equations of step 2 are applied to plot the 
Figures (see Appendix D). 
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Figure 6: e = 1, a = 1, A = 0.4 and k = 1.3 
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Figure 7: e — 1.6, a — 1.3, A = 0.9 and k — 
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Figure 8: e = 0.8, a = 0.2, A = 0.6 and k = 
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Figure 10: e = 0.6, a = 1.2, A = 1.2 and k = 1.7 
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Figures 11 15 in these profiles of the normalized concentrations of the sub- 
strate u, enzyme-substrate complex v, enzyme E and product w correspond 
the case 1-case 5, respectively. The equations of step 3 are applied to plot the 
Figures (see Appendix E). 




Figure 11: e = 1.001, a = 1.001, A = 0.4 and k = 1.3 
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Figure 12: e = 1.6, a = 1.3, A = 0.9 and k = 1.7 
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Figure 15: e = 0.6, a = 1.2, A = 1.2 and k = 1.7 



Figures |16}j20| in these profiles of the normalized concentrations of the sub- 
strate u, enzyme-substrate complex v, enzyme E and product w correspond 



the case 1-case 5, respectively. The equations ( 31 )-( 35 ] are applied to plot the 
Figures (see Appendix A). 




Figure 16: e = 1, a = 1, A = 0.4 and k = 1.3 
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Now a days, there are many attempts to find the approximate solutions in the 
filed of biochemical reactions, especially for those problems which are difficult 
to find the exact solution of them and to describe the behavior of their reactions 
and parameters. These attempts will probably continues to approach the best 
way to link the theoretical studying such as kinetic equations with experimental 
results. The system of reactions |l]) is one of the case study in the past which 
was developed in different ways. Then, the system of reactions ^ is more 
general than ([I]) in terms of reactions and it is not easy to find exact solution 
of them. So, our results in this section based on system ^ which is described 
in system of equations ( 18 )-( 19 ) . We have applied two different methods for 
equations (18)-(19l to find the analytical approximate solutions. The first one 
(HPM) which was used by many researchers for system |T]) [231 [22 [2U E] , and 
the second method (SIM) is our simple approach to find anylatical approximate 
solution for the same problem. The main aim of this discussion is to emphasis 
that there is a wide range data results in our work to show the goodness of the 
approximate solution of equations (|18h— ( 19 ). Another aim of this section is that 



to recognize the best iteration of (SIM) compared to (HPM). 
There is a variety of data results that is why our approach (SIM) is quite close 
from the classical method (HPM) to find the approximate solution in terms 
of the value of concentrations u,v,E and w. First of all, the second iteration 
in the (SIM) has many significant similarities with (HPM) and some of them 
have a grate result of our work, for instance, the Figures |6fT0] show that the 
value of concentration of substrate u slightly decease from its initial value the 
concentration (u(0) = 1) and there is a few changes in the value of concentra- 
tion of enzyme-substrate complex v. Generally, they reach to some constant 
values after r > 4. Also, in the Figures [TBpO] appear that the rate of these 
components are slightly similar with the corresponding Figures |6][T0] and they 
reach to some constant values in each cases. Another example is that the value 
of concentration of enzyme E in both set of figures more or less has the same 
value of concentrations, especially in case 1, case 2 and case 5. The similarities 
are occurred in the same cases for value of concentration product w. 
Another interesting point is that in our technique Figure[l3]in step 3 the value of 
concentration v reaches to max point when < r < 2, also in the same interval 
of time the value of concentration v reaches to max as in Figure 18 as well. We 
can also realize that the value of enzyme in both figures ends up the min value 
when < r < 2. In addition, the Figures [TTfTS] and Figures [TBpO] show that 
there is a gradual decrease in the rate of substrate u between < t < 2 and it 
reaches to level off after r > 4. Whereas, the rate of product w slightly increase 
of both figures between < r < 2 and it is likely to remain stable after r > 4. 
Generally speaking, it is more reliable to see hole similarities in terms of their 
values by comparison the figures. As a result, it could be argued that the 
iterations in the simple approach (SIM), especially the last two iterations are 
more similar to find the approximate solution for system equations (18|-(19| in 
compared to (HPM). 

On the other hand, there is some differences between our simple technique 
(SIM) and classical technique (HPM) in in term of approximate solution, for 
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example, the Figures [T][5] show that the value of concentration of substrate u 
slightly decease from its initial value the concentration (u(0) = 1) and there 
is a few changes in the value of concentration of enzyme-substrate complex v. 
Generally, they become zero after r > 5. While, in the Figures [T6p0| appear 
that the rate of the components are not become zero value but they reach to 
some constant values in each case. Basically, It could be pointed out that the 
differences between of them are small and it can be ignored. 
Overall, it seems that the data results of our study are clear examples to support 
our simple technique compared to the (HPM) and there is so many similarities 
between of them in terms of their rates, for example, the similarities are quite 
clear in the Figures [6]|T0| and the Figures [TBpOl Although, it could be say that 
there is some different results in their values, but it is occurred as a tiny value. 
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7 Findings 



From the previous section, we have applied two methods to estimate the ap- 
proximate solution to our case study, and it could be argued that the methods 



are quite appropriate to describe the parameters in our problem equations (181 



(19). Although, it is not easy to compare all results in the previous section to 
estimate the appropriate iteration in (SIM) in one side and (HPM) in another 
side. Because of this, from the Tables fT|3] we have used the mean of the second 



norm equation (47) to find the total differences between (HPM) and each itera- 



tion of (SIM). Then, in the Figure 21 we point out that the rate of convergence 
between the (SIM) and (HPM). This formula helps us to make sure about our 
results in this paper. Therefore, we use the following equation to find the results 
in the Tables [HI 

^/W~ = V n (47) 

Where / and g are the value of concentration of substances u, v, E and w for 
(SIM) and (HPM), respectively, and N is the number of timescale iterations. 

The second iteration of (SIM) is more similar to (HPM) regarding our find- 
ings in Figure [2l] (see H-S2 values). On the other hand, the differences between 
our approach (SIM) and (HPM) are more occurred in the first iteration than 
other iterations (see H-Sl values). It could be said that this iteration is not 
quite appropriate in this case study. It may be caused by putting the non-linear 
part in this iteration as zero value (see step 1). Whereas, there is a small value 
of average norm between second iteration and (HPM). For instance, the average 
value of norm concentration of E is just a small value and equal to 0.02 (sec 
H-S2 for NE value). It means, the second iteration method is more a appro- 
priate iteration one in this work to approach the approximate solution. Also, 
this small value in terms of the mean of norm is occurred in third iteration as 
well (see H-S3 values). Generally speaking, the second iteration method of our 
work is the best iteration one to get the convergence in terms of the solution in 
comparison with the classical method (HPM). 
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Norm Con. 


Case 1 


Case 2 


Case 3 


Case 4 


Case 5 


Average 1 


UU1 


0.26747 


0.20151 


0.17482 


0.16551 


0.18084 


0.19803 


VV1 


0.26918 


0.20313 


0.17535 


0.16627 


0.18385 


0.199556 


EE1 


0.29684 


0.24889 


0.14721 


0.18732 


0.31446 


0.238944 


WW1 


0.5633 


0.37454 


0.36304 


0.29118 


0.614 


0.441212 



Table 1: shows that the average number of second norm between the first itera- 
tion of (SIM) and (HPM) for five different cases in terms of different parameters 
by using the formula of equation ( 47 1 . It can be easily calculated the results of 
the Table [I] by using Matlab program in (Appendix B and C). 



Nor. Con. 


Case 1 


Case 2 


Case 3 


Case 4 


Case 5 


Aver. 2 


UU2 


0.03611 


0.02504 


0.02326 


0.14113 


0.05728 


0.056564 


VV2 


0.03611 


0.02507 


0.02333 


0.14114 


0.05728 


0.056586 


EE2 


2.204 E-17 


0.01363 


0.02554 


0.06548 


0.01727 


0.024384 


WW2 


0.03611 


0.02008 


0.05259 


0.16332 


0.08771 


0.071962 



Table 2: shows that the average number of second norm between the second 
iteration of (SIM) and (HPM) for five different cases in terms of different pa- 



rameters by using the formula of equation ( 47 ) . It can be easily calculated the 
results of the Table [2] by using Matlab program in (Appendix B and D). 



8 Conclusion 

Simple iteration method (SIM) and Homotopy perturbation method (HPM) are 
used to approximate analytic solutions to the non-linear differential equations. 
A simple, straight forward methods are derived for estimating the concentrations 
of substrate u, product w, enzyme-substrate complex v and enzyme E. The 
dimensionless technique can be simply generalized to all kinds of system of 
non-linear differential equations with numerous complex boundary conditions 
in enzyme-substrate reaction diffusion processes |21j . 

Homotopy perturbation method (HPM) has been used for simple enzyme 
reaction ([I]) and it gives a good approximate solution [TTJ [5T] . We have applied 



Norm Con. 


Case 1 


Case 2 


Case 3 


Case 4 


Case 5 


Average 3 


UU3 


0.08435 


0.04113 


0.04704 


0.13439 


0.03628 


0.068638 


VV3 


0.08439 


0.04113 


0.0471 


0.13448 


0.03633 


0.068686 


EE3 


0.02625 


0.00942 


0.02733 


0.0672 


0.03269 


0.03257 


WW3 


0.10839 


0.02597 


0.08244 


0.15919 


0.06823 


0.088844 



Table 3: shows that the average number of second norm between the third itera- 
tion of (SIM) and (HPM) for five different cases in terms of different parameters 
by using the formula of equation ( 47 1 . It can be easily calculated the results of 
the Table [| by using Matlab program in (Appendix B and E). 
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■ H-Sl 

■ H-S2 

H-S3 



u y e w 



Figure 21: shows that the average value of second norm convergence of the value 
of concentrations u,v,E and w between the Homotopy perturbation method 
(HPM) and correspond the iterations of simple iteration method (SIM). The 
blue column (H-Sl), brown column (H-S2) and green column (H-S3) are de- 
scribed the second norm differences between (HPM) and the first, second and 
third iterations, respectively. We use Excel to plot this figure based on the 
results in Tables Q]|3] 

. The symbols u, V, e and w in this figure just represent the second norm 
differences between the two methods which are applied in this project. 



this method for our chemical reaction ^ and we have obtained the approximate 
solution such as other people did but for different project. Furthermore, we have 
derived a simple approach which is called Simple iteration method (SIM), this 
method is based on some iterations (steps) in the second step the approximate 
solution is quite close to the classical method (HPM), because the value of all 
concentrations in each case are similar with the classical one (see Figures 6fl0 
and Figures [16J20 1. 

We have also used the idea of second norm to point out the best iteration 
for our approximation. So, it is clear that the second iteration method has quite 
similar with (HPM) to obtain the solution. Consequently, the Figure 21 includes 
a golden point that why our simple technique is important and has a good 



similarities with (HPM), particularly the second iteration method (see Figure 21 



H-S2 values). Basically, it could be applied this simple technique (SIM) for some 
other complex chemical reactions to find appropriate solution of them and to 
describe the behavior of their parameters much more easier than methods which 
were used in many recent works. For example, it will be applied to many open 
path ways of biochemical reactions [3 . In addition, we will highly recommend to 
apply the simple approach (SIM) to describe approximate solution of complex 
enzyme reactions which are not easy to find out their solutions |14) , and a future 
study investigating for the reaction mechanism of competitive inhibition and the 
reaction scheme of allosteric inhibition would be very interesting 0]. 
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Appendix A This appendix consists of the solution of equations ( 18 1 and 
(19)) by using Homotopy perturbation method. Furthermore, this method is 
used to drive equations (31 1 and (32) from equations (18) and (19), let a = 
e(k — A), b = 1 — a and c = k + a + as, 



lA .,du . ,du 
(1 - — I" £u \ + Q-l~, — \- £u — av - 
dr dr 



0, 



(48) 



dv du 

(1 ~~ — I" cv \ + bu + cv + buv — aev 2 — a] = 0, (49) 

dr dr 



with initial conditions, 



u(0) = 1,?;(0) = 



(50) 



Thus, by Homotopy perturbation method [TT] [5T1 [52] , the approximate solution 
of equations ( 48 ) and ( 49 1 are: 



u = Uq + qui + q u 2 + 
v = v + qvi + q 2 v 2 + ■ 



(51) 
(52) 
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Substituting the equations (51) and (52 1 in the equations (48) and (49) and 
comparing the coefficient of like power q, we can get the following system of 
ordinary differential equations: 



and, 



q v : — + eu = 
dr 

q : — h eu\ — avo — euqVq = 

dr 

2 dU2 

q : — h eu2 — avi — euqVi — euiVq = 

dr 



o dv o . n 
q ■ -r- + cvq = 
dr 

q : — h cvi — oito + ouo u o ~ ae^o — a = 0, 

(XT 

q 2 : — — + CV2 — bu\ + biiQVi + bu\Vo — 2ai?vo w i = 0, 
dr 



(53) 

(54) 
(55) 

(56) 
(57) 
(58) 



by solving the equations ( 53)-( 58), with using initial conditions equation (50), 
we obtain the following solutions: 



u 2 (t) 



and, 



ab ae 

c — e c 



te- 



as — cbe — cae 



uo(r) = e~ ST 

u 1 (r) = e- £T 

abc — acts + aac 
c(s — c) 2 
act b 



(59) 
(60) 



(-e-c)r 



-2er 



c 2 (e — c) ce 
acts — abc — aca act 



-1 + 



e — c 

b —ae 2 + cbe + cae 



c(e — c) 5 



ce c — e 



c 2 (e — c) 



(61) 



voir) - 



v 2 (t) = 



v 1 (t) = e 

c — e 



ba 



c — e c(c — e) 
aeb — cb 2 — cba 
ce(e — c) 
b 



be — ae + c 

H 

c(e — c) 

b 2 

(c-e)(c-2e) ( 



-e 
b 2 



(-e-c)t 



cb 2 — bae + cba 



ab 



e — c (e — c)(c — 2e) ce(e — c) c(e — c) 



(62) 
(63) 



(64) 
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According to Homotopy perturbation method (HPM), we can easily find that 



u(t) — lim u(t) = Uq + U\ + u 2 + 



(65) 



and, 



v(t) 



lim v(t) — v + vi + v 2 + 

q->l 



(66) 



( 66 ) , we obtain the approximate solution for system of non-linear ODE equations 



We put equations ( 59 )-( 61 1 in equation ( 65 1 and equations ( 62 )-( 64 ) in equation 



( 18 )-( 19 1 which is described in equations (31) and (32) 



Appendix B Let kl — e, k2 — A, fc3 = a and t — r. then we use the following 
Matlab programming to plot the functions of equations (31)-(35), and we have 
shown the the affect of the parameters on the chemical substances u, v, E and 
w in Figures [TBpOl 
t = 0; 

for i = 1:101 

kl = 1 ; k2 = 1.2; k3 = 0.9; k = 1.3; 

a = kl*(k-k2); b = l-k3; c = k+k3*kl+k3; 

u = 2*exp(-kl*t)+((a*b)/(c-kl) + (k3*kl)/c)*t*exp(-kl*t) 

+ (a*b*c-a*k3*kl+a*k3*c)/(c*(fcl - c) 2 )*exp(-c*t) 

+ b/(kl-c)* exp(-2*kl*t)+(k3*fcl 2 -c*b*kl-k3*c*kl)/(c 2 *(kl-c)) 

*exp((-kl-c)*t)+ (a*k3)/(c*kl) 

+ (-l+(a*kl*k3-a*b*c -a*k3*c)/(c*(fcl - c) 2 )-(a*k3)/(c*kl) 

+ b/(c-kl) + (a*b*kl-k3*A:l 2 +k3*c ,,! kl)/(c 2 *(kl-c)))*exp(-kl !,t t); 

v = b/(c-kl)*exp(-kl*t)+((c*b-k3*kl+k3*c)/(c*(kl-c))) !,c exp(-c*t)+k3/c 

+ (b/(c-kl) + (k3*b)/(c*(c-kl)))*exp(-kl*t) 

+ (fe 2 /((c-kl)*(c-2*kl)))*exp(-2*kl*t) 

+(k3*kl*b-c*6 2 -k3*c !,t b)/(c*kl 5,! (kl-c)) 5,! exp((-c-kl)*t) 

+ (b/(kl-c)+6 2 /((kl-c)*(c-2*kl)) + (c*6 2 -k3*kl*b 

+k3*c*b)/(c*kl*(kl-c)) + (k3*b)/(c*(kl-c)))*exp(-c 5,! t); 

E = 1-v; w = 1/kl -u/kl -v; A(i) = u; B(i) = v ; C(i) = E; D(i) = w; T(i) = t; 

t = t+0.1; 

end 

plot(T,A,'r',T,B,'k.',T,C,'b.',T,D,'g') 
ylabel('Concentration of u, v, E and w') 
xlabel('Dimensionless Time(t)') 
axis square 



Appendix C Let kl = e, k2 — A, kS — a and t — r, then we use the following 
Matlab programming to plot the functions of step 1, and we have shown the 
the affect of the parameters on the chemical substances u, v, E and w in Figures 

he 

t = 0; 

for i = 1:101; 

kl = 1; k2 = 1.2; k3 = 0; k = 1.3; a=kl*(k-k2); b=l-k3; c=k+k3+k3*kl; 
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pi = (-kl-c-sqrt((M+c) 2 -4*(c*kl-a*b)))/2 ; p2 = (-kl-c+sqrt((fcl+c) 2 -4*(c*kl- 
a*b)))/2 ; 

d2 = ( P 2+kl)/(p2-pl); d3 = (pl+kl)/(pl- P 2); 
dl = ((pl+kl)*(p2+kl))/(a*(p2-pl)); 

ul = d2*exp(pl*t)+d3*exp(p2*t); vl = dl*exp(pl*t)-dl*exp(p2*t); 

el = 1-vl; wl = 1/kl -ul/kl -vl; 

A(i) = ul; B(i) = vl; C(i) = el; D(i) = wl; T(i) = t; 

t = t+0.1; 

end 

plot(T,A,'r',T,B,'k.',T,C,'b.',T,D,'g') 
ylabel('Concentration of wl, ul, SI and wV) 
xlabel('Dimensionless Time(t)') 
axis square 



Appendix D Let kl = e, k2 = A, k'i = a and t = r, then we use the following 
Matlab programming to plot the functions of step 2, and we have shown the 
the affect of the parameters on the chemical substances u,v,E and w in Figures 

mm 

t = 0; 

for i = 1:101; 

kl = 1; k2 = 1.2; k3 = 0.2; k = 1.3; a = kl*(k-k2); b = l-k3; c = k+k3+k3*kl; 

pi = (-kl-c-sqrt((fcl + c) 2 -4*(c*kl-a*b)))/2 ; 

p2 = (-kl-c+sqrt((fcl + c) 2 -4*(c*kl-a*b)))/2 ; 

d2 = ( P 2+kl)/(p2-pl); d3 = (pl+kl)/(pl- P 2); 

dl = ((pl+kl)*(p2+kl))/(a*(p2-pl)); 

d4 = klWdl; d5 = (-kl*d2*dl)+(kl*d3*dl); 

d6 = -kl*d3*dl; d7 = (-b*d2*dl) + (k3*kl*dl 2 ; 

d8 = (b*d2*dl)-(b*d3*dl)-(2*k3*kl*dl 2 ); 

d9 = (b*d3*dl)+(k3*kl*dl 2 ) ; dlO = ( P 2+kl)/(a*(p2-pl)); 

dll = -l/( P 2-pl); dl2 = (pl+kl)/(a*(pl-p2)); 

dl3 = l/(p2-pl); dl4 = (d4*dl0+dll*d7)/pl ; 

dl5 = (dl0*d5+dll*d8)/p2 ; dl6 = (dl0*d6+dll*d9)/(2*p2-pl); 

dl7 = (-k3*dll)/pl ; dl8 = (dl2*d4+dl3*d7)/(2*pl- P 2); 

dl9 = (dl2*d5+dl3*d8)/pl ; d20 = (dl2*d6+dl3*d9)/p2 ; 

621 = (-k3*dl3)/p2 ; d22 = a*dl4+a*dl8; d23 = a*dl5+a*dl9; 

d24 = a*dl6+a*d20; d25 = a*dl7+a*d21; d26 = ( P l+kl)*dl4+(p2+kl)*dl8; 

d27 = (pl+kl)*dl5+(p2+kl)*dl9 ; d28 = ( P l+kl)*dl6+(p2+kl)*d20 ; 

d29 = ( P l+kl)*dl7+(p2+kl)*d21 ; 

c3 = (l-d22-d23-d24-d25)/a -(a*(d26+d27+d28+d29) 

+ ( P l+kl)*(l-d22-d23-d24-d25))/(a*(pl-p2)); 

c4 = (a*(d26+d27+d28+d29) 

+ ( P l+kl)*(l-d22-d23-d24-d25))/(a*(pl-p2)); 

u2 = c3*a*exp(pl*t)+c4*a*exp(p2*t)+d22*exp(2*pl*t) 

+d23*exp((pl+p2)*t)+d24*exp(2*p2*t)+d25; 

v2 = c3*(pl+kl)*exp(pl*t)+c4*(p2+kl)*exp(p2*t) 
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+d26*exp(2*pl*t)+d27*exp((pl+p2)*t)+d28*exp(2*p2*t)+d29; 

c2 = l-v2; w2 = 1/kl -u2/kl -v2; 

A(i) = u2; B(i) = v2; C(i) = c2; D(i) = w2; T(i) = t; 

t = t+0.1; 

end 

plot(T,A,V,T,B,'k.',T,C,'b.',T,D,'g') 
ylabel('Concentration of u2,v2,E2 and w2') 
xlabel('Dimensionless Time(i)') 
axis square 

Appendix E Let fcl = e, k2 = A, fc3 = a and t = t, then we use the following 
Matlab programming to plot the functions of step 3, and we have shown the 
the affect of the parameters on the chemical substances u, v, E and w in Figure 

EEi 

t = 0; 

for i = 1:101; 

kl = 1; k2 = 1.2; k3 = 0.2; k = 1.4; a = kl*(k-k2); 

b = l-k3; c=k+k3+k3*kl; 

pi = (-kl-c-sqrt((fcl + c) 2 -4*(c*kl-a*b)))/2; 

p2 = (-kl-c+sqrt((fel + c) 2 -4*(c*kl-a*b)))/2 ; 

d2 - ( P 2+kl)/(p2-pl); d3 - (pl+kl)/(pl- P 2); 

dl = ((pl+kl)*(p2+kl))/(a*(p2-pl)); d4 = kl*d2*dl; 

d5 = (-kl*d2*dl)+(kl*d3*dl); 

d6 = -kl*d3*dl; d7 = (-b*d2*dl) + (k3*kl*dl 2 ); 

d8 = (b*d2*dl)-(b*d3*dl)-(2*k3*kl*dl 2 );d9 = (b*d3*dl)+(k3*kl*dl 2 ); 

dlO = (p2+kl)/(a*(p2-pl)); dll=-l/(p2-pl); dl2 = (pl+kl)/(a*(pl-p2)); 

dl3 = l/( P 2-pl); dl4=(d4*dl0+dll*d7)/pl ;dl5 = (dl0*d5+dll*d8)/p2 ; 

dl6 = (dl0*d6+dll*d9)/(2*p2-pi); dl7 = (-k3*dll)/pl; 

dl8 = (dl2*d4+dl3*d7)/(2*pl-p2); dl9 = (dl2*d5+dl3*d8)/pl ; 

d20 = (dl2*d6+dl3*d9)/p2 ; d21 = (-k3*dl3)/p2 ; d22 = a*dl4+a*dl8; 

d23 = a*dl5+a*dl9; d24 = a*dl6+a*d20; d25 = a*dl7+a*d21; 

d26 = ( P l+kl)*dl4+(p2+kl)*dl8; d27 = (pl+kl)*dl5+(p2+kl)*dl9 ; 

d28 = ( P l+kl)*dl6+(p2+kl)*d20 ; d29 = ( P l+kl)*dl7+(p2+kl)*d21 ; 

c3 = (l-d22-d23-d24-d25)/a-(a*(d26+d27+d28+d29) 

+ ( P l+kl)*(l-d22-d23-d24-d25))/(a*(pl-p2)); 

c4 = ( a *(d26+d27+d28+d29)+(pl+kl)*(l-d22-d23-d24-d25))/(a*(pl-p2)); 
d30 = kl*c3 2 *a*(pl+kl)+kl*d22*d29+kl*d26*d25; 

d31 = kl*c3*c4*a*(p2+kl)+kl*c4*c3*a*(pl+kl)+kl*d23*d29+kl*d25 !,t d27; 
d32 = kl*c3*a*d26+kl*c3*(pl+kl)*d22; 

d33 = kl*c3*a*d27+kl*c4*a*d26+kl*c4*(p2+kl)*d22+kl*c3*(pl+kl)*d23; 

d34 = kl*c3*a*d28+kl*c4*a*d27+kl !,! c4*(p2+kl)*d23+kl*c3*(pl+kl) !,! d24; 

d35 = kl*c3*a*d29+kl*c3*(pl+kl)*d25; 

d36 = kl*c4 2 *a*(p2+kl)+kl*d24*d29+kl*d25*d28; 

d37 = kl*c4*a*d28+kl*c4*(p2+kl*d24; 

d38 = kl*c4*a*d29+kl*c4*(p2+kl)*d2; 
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d39 = kl*d22*d26; d40 = kl*d22*d27+kl*d23*d26; 

d41 = kl*d22*d28+kl*d23*d27+kl*d24*d26; d42 = kl*d23*d28+kl*d24*d27; 
d43 = kl*d24*d28; d44 = kl*d25*d29; 

d45 = -b*d30/kl+k3*kl*c3 2 *(pl + fcl) 2 +k3*kl*d26*d29+k3*kl*d29*d26; 

d46 = -b*d31/kl+k3*kl*c3*c4*(pl+kl)*(p2+kl) 

+ k3%l*c3*c4*(pl+kl)*(p2+kl)+k3*kl*d27*d29+k3*kl*d27*d29; 

d47 = -b*d32/kl+k3*kl*c3*d26*(pl+kl)+k3*kl*c3*d26*(pl+kl); 

d48 = -b*d33/kl+k3*klWd27*(pl+kl)+k3*kl*c4*d26*(p2+kl); 

d49 = -b*d34/kl+k3*kl*c3*d28*(pl+kl)+k3*kl*c4*d27*(p2+kl) 

+ k3%l*c4*d27*(p2+kl)+k3*kl*c3*d28*(pl*kl); 

d50 = -b*d35/kl+k3*kl*c3*d29*(pl+kl)+k3*kl*c3*d29*(pl*kl); 

d51 = -b*d36/kl+k3*kl*c4 2 *(p2+kl)+k3*kl*d28*d29+k3*kl*d28*d29; 

d52 = -b*d37/kl+k3%l*c4*d28*(p2+kl)+k3*kl*c4*d28*(p2+kl); 

d53 = -b*d38/kl+k3*kl*c4*d29*(p2+kl)+k3*kl*c4*d29*(p2+kl); 

d54 = -b*d39/kl+k3*kl*(d26) 2 ; 

d55 = -b*d40/kl+k3*kl*d26*d27+k3*kl*d26*d27; 

d56 = -b*d41/kl+k3*kl*d26*d28+k3*kl*(d27) 2 +k3*kl*d26*d28; 

d57 = -b*d42/kl+k3*kl*d27*d28+k3*kl*d27*d28; 

d58 = -b*d43/kl+k3*kl*(d28) 2 ; d59 = -b*d44/kl+(d29) 2 +k3; 

d60 = (dl0*d30+dll*d45)/(pl); d61 = (dl0*d31+dll*d46)/(p2); 

d62 = (dl0*d32+dll*d47)/(2*pl); d63 = (dl0*d33+dll*d48)/(pl+p2); 

d64 = (dl0*d34+dll*d49)/(2*p2); d65 = (dl0*d35+dll*d50); 

d66 = (dl0*d36+dll*d51)/(-pl+2*p2); 

d67 = (dl0*d37+dll*d52)/(-pl+3*p2); 

d68 = (dl0*d38+dll*d53)/(-pl+p2); d69 = (dl0*d39+dll*d54)/(3*pl); 
d70 = (dl0*d40+dll*d55)/(2*pl+p2); d71 = (dl0*d41+dll*d56)/(pl+2*p2); 
d72 = (dl0*d42+dll*d57)/(3*p2); d73 = (dl0*d43+dll*d58)/(-pl+4*p2); 
d74 = (dl0*d44+dll*d59)/(-pl); d75 = (dl2*d30+dl3*d45)/(2*pl-p2); 
d76 = (dl2*d31+dl3*d46)/(pl); d77 = (dl2*d32+dl3*d47)/(3*pl-p2); 
d78 = (dl2*d33+dl3*d48)/(2*pl); d79 = (dl2*d34+dl3*d49)/(pl+p2); 
d80 = (dl2*d35+dl3*d50)/(pl-p2); d81 = (dl2*d36+dl3*d51)/(p2); 
d82 = (dl2*d37+dl3*d52)/(2*p2); d83 = (dl2*d38+dl3*d53); 
d84 = (dl2*d39+dl3*d54)/(4*pl-p2); d85 = (dl2*d40+dl3*d55)/(3*pl); 
d86 = (dl2*d41+dl3*d56)/(2 !,! pl+p2); d87 = (dl2*d42+dl3*d57)/(pl+2*p2); 
d88 = (dl2*d43+dl3*d58)/(3*p2); d89 = (dl2*d44+dl3*d59)/(-p2); 
d90 = a*d60+a*d75;d91=a*d61+a*d76; 
d92 = a*d62+a*d77; d93=a*d63+a*d78; 

d94 = a*d64+a*d79; d95 = a*d65; d96 = a*d80; d97 = a*d66+a*d81; 
d98 = a*d67+a*d82; d99 = a*d68; dlOO = a*d83; dlOl = a*d69+a*d84; 
dl02 = a*d70+a*d85; dl03 = a*d71+a*d86; dl04 = a*d72+a*d87; 
dl05 = a*d73+a*d88; dl06 = a*d74+a*d89; hi = pl+kl; h2 = p2+kl; 
dl07 = hl*d60+h2*d75; dl08 = hl*d61+h2*d76; dl09 = hl*d62+h2*d77; 
dllO = hl*d63+h2*d78; dill = hl*d64+h2*d79; 
dll2 = hl*d65; dll3 = h2*d80; 
dll4 = hl*d66+h2*d81; dll5 = hl*d67+h2*d82; 
dll6 = hl*d68; dll7 = h2*d83; 
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dll8 = hl*d69+h2*d84; dll9 = hl*d70+h2*d85; dl20 = hl*d71+h2*d86; 

dl21 = hl*d72+h2*d87; dl22 = hl*d73+h2*d88; dl23 = hl*d74+h2*d89; 

M = d90+d91+d92+d93+d94+d96+d97+d98+d99 

+dl01+dl02+dl03+dl04+dl05+dl06; 

N = dl07+dl08+dl09+dll0+dlll+dll3+dll4+dll5 

+dll6+dll8+dll9+dl20+dl21+dl22+dl23; 

c5 = ((l-M)*h2+a*N)/(a*(h2-hl)); c6=-N/h2-(hl/h2)*c5; 

u3 = c5%*exp(pl*t)+c6*a*cxp(p2*t)+d90*exp((2*pl)*t) 
+d91*oxp((pl+p2)*t)+d92*exp((3*pl)*t) 

+d93*exp((2*pl+p2)*t)+d94*exp((pl+2*p2)*t)+d95*t*exp((pl)*t) 

+d96*exp((pl)*t)+d97*exp((2*p2)*t)+d98*exp((3*p2)*t) 

+d99*exp((2*p2)*t)+dl00*t*cxp((2*p2)*t)+dl01*cxp((4*pl)*t) 

+dl02*exp((3*pl+p2)*t)+dl03*exp((2*pl+2*p2)*t) 

+dl04*exp((pl+3*p2)*t)+dl05*cxp((4*p2)*t)+dl06; 

v3 = c5%l*exp(pl*t)+c6*h2*exp(p2*t)+dl07*exp((2*pl)*t) 

+dl08*exp((pl+p2)*t)+dl09*exp((3*pl)*t) 

+dll0*exp((2*pl+p2)*t)+dlll s,! exp((pl+2*p2)*t)+dll2*t*cxp((pl)*t) 

+dll3*exp((pl)*t)+dll4*exp((2*p2)*t)+dll5*exp((3*p2)*t) 

+dll6*exp((2*p2)*t)+dll7*t*exp((2*p2)*t) 

+dll8*exp((4*pl)*t)+dll9*cxp((3*pl+p2)*t)+dl20*cxp((2*pl+2*p2) !,c t) 

+dl21*cxp((pl+3*p2)*t)+dl22*exp((4*p2)*t)+dl23; 

c3 = l-v3; w3 = 1/kl -u3/kl -v3; 

A(i) = u3; B(i) = v3; C(i) = c3; D(i) = w3; T(i) = t; 

t = t+0.1; 

end 

plot(T,A,'r',T,B,'k.\T,C,'b.',T,D,'g') 
ylabel('Concentration of u3,v3,E3 and u>3') 
xlabel('Dimcnsionless Time(i)') 
axis square 
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